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1.  INTRODUCTION 


1  2 

Preceding  reports  ’  by  the  authors  describe  an  algorithm  for  computing  derivatives  of 
a  function  of  a  real  variable  from  its  tabular  data.  The  algorithm  is  based  on  two  sets  of 
assumptions.  The  first  set  specifies  a  class  of  functions  such  that  the  tabular  data  are 
approximate  values  of  one  of  these  functions.  The  second  set  of  assumptions  determines 
a  finite  dimensional  subspace  Ak  of  the  selected  functions  from  which  an  approximation  of 
the  tabular  data  is  selected.  The  essential  difference  between  the  current  algorithm  (i.e., 
the  algorithm  described  in  this  report)  and  any  of  the  standard  procedures  such  as  moving 
polynomial  arc,  linear  regression,  spline  approximation,  and  others  is  due  to  this  second 
set  of  assumptions.  The  standard  methods  assume  directly  or  indirectly  that  an  element  of 
a  preselected  finite  dimensional  space,  say  Sk ,  yields  an  appropriate  approximation  and  that 
Sk  is  independent  of  the  data  at  hand.  The  "best"  approximation  in  Sk  is  obtained  by 
applying  a  preselected  criterion  characteristic  of  the  algorithm.  In  contrast  to  this  our  pro¬ 
cedure  examines  the  data  and  selects  a  linear  subspace  Ak  of  an  assumed  infinite  dimen¬ 
sional  algebra  A  that  is  most  appropriate  for  the  available  data  and  then  determines  a  fam¬ 
ily  of  approximations  in  the  subspace  Ak  dependent  on  the  data  and  also  the  relative  accu¬ 
racy  of  these  approximations  expressed  by  weights  which  in  their  turn  depend  on  the  data. 
Final  results  are  weighted  averages  of  individual  approximations. 

Usually  a  subspace  Sk  and  its  approximating  element  are  selected  by  a  heuristic  cri¬ 
terion  such  as  minimum  root  mean  square  error,  degree  of  smoothness,  and  others.  In 
most  cases  this  criterion  represents  a  compromise  between  simplicity  (smoothness,  dimen¬ 
sion  of  Sk )  of  an  approximating  function  and  its  faithfulness  to  the  data.  Similarly,  our 
procedure  selects  approximating  elements  and  their  weights  by  a  compromise  between  the 
accuracy  of  approximation  and  the  robustness  of  the  approximating  model.  The  criterion 
for  this  compromise  is  purely  heuristic.  In  this  report  we  examine  such  a  criterion,  which 
is  a  modification  of  one  used  in  the  previous  reports.  ’ 

2 

We  note  that  approximating  functions  selected  by  the  algorithm  yield  no  values  of 
derivatives  in  the  initial  and  final  segments  of  the  data,  where  the  lengths  of  these  seg¬ 
ments  are  dependent  on  the  data.  The  number  of  points  lost  in  this  way  depends  on  the 
dimension  k  of  the  subspace  Ak  and  on  the  selected  multiple  q  of  the  data  step  size.  Here 
we  describe  a  procedure  for  computing  derivatives  at  the  beginning  and  the  end  of  the 
data  sequence.  These  are  obtained  with  the  aid  of  digital  linear  filters  dependent  on  the 
data. 


1  C.  Masaitis  and  G.  Francis,  " Numerical  Differentiation  of  Noisy  Data,"  ARBRL-MR-03126,  Aug  81,  AD 
A104631. 

2  C.  Masaitis  and  G.  Francis,  ” Computation  qf  Derivatives  from  Tabular  Data,"  ARBRL-MR-03188,  Jul 
82,  AD  B066124L. 
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2.  ALGORITHM 


Details  of  the  algorithm  for  computing  derivatives  are  described  in  the  report 1  which  is 
denoted  by  ND  in  the  references  below.  A  summary  of  these  details  follows. 

Various  formulas  of  this  algorithm  are  derived  from  four  assumptions: 

a.  Tabular  data  x(n),n  —  1,2 ,..JV  are  obtained  by  measuring  a  function  yit )  at 
equally  spaced  points  with  a  step  size  A.  The  measuring  error  e„  is  a  white  noise  with  zero 
mean  and  unknown  variance  <x2  independent  of  time,  i.e.  x  in  )  —  y  in  A)  +  en . 

b.  The  function  yit )  belongs  to  a  family  of  functions  whose  derivatives  on  the 
interval  [0,T]  can  be  represented  by  a  linear  combination  of  values  of  y(f ),  i.e. 

y'(f)  —  £  at(y)y[t  +(v  —  /)  A  ]  (2.1) 

/  -  1 

for  any  integer  v  satisfying  —t+k  A^v  A  —  t.  The  coefficients  a,(v)  depend  on 
the  function y(f)  and  the  location  parameter  v  but  not  on  t  . 

As  shown  in  Section  3  of  ND  assumption  ib  )  implies  that  y  it)  belongs  to  the  algebra 
A  generated  by  algebraic,  trigonometric,  and  exponential  polynomials  defined  on  the  inter¬ 
val  [0,7"  ]. 

According  to  Section  1  of  ND  an  assumption  that  yit)  €  A  is  equivalent  to  an  assump¬ 
tion  that  for  every  A  >0  yit)  satisfies  the  difference  equation: 

P  A  (5  Jy  (/)-0,  (2.2) 

where  P  &  (X)  —  ft  ~  +  lfmd  y  (O—y  (t  —  A).  Hence,  an  approximating 

7-1 

function  can  be  obtained  by  constructing  an  autoregressive  model  such  as  (2.2). 

We  write  x  in,p,q)  for  xip+qn),  for  any  positive  integer  q  and 

N  —  p 

p  =0,  1,  2,..., q  -  l,#i  =  0,1, ...JVp,  where  Np  —  [ — - — 1,  i.e. 

x  {n,p,q)  =yir  A^e,.,  (2.3) 

where  r  —  p  +  qn .  By  substituting  this  in  (2.2)  we  get 

x  (n,p,q)  — er  -  £  [x  (n-y ,/>,<?)- €r_;l.  (2.4) 

7  —  1  1  J 

According  to  Section  2  of  ND  the  minimum  variance  estimates  of  the  fly's  are  obtained  by 
the  following  iterative  procedure. 
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Let  M  (k,q)  be  the  matrix  of  the  normal  equations  of  the  overdetermined  system 


k 

£  aj  x  q)-=x  (n,p,  q),  (2.5) 

7-  1 


n  —  k  +  1,...,  N„,p  —  0,  1  ,...,?  —  1.  Let  X  (k,q)  be  the  right-hand  side  vector  of  these 

r  ^  | 

normal  equations  and  N  —  £  {Np  —  k  ).  Let 

p  -0 


1  *£' 


ar  u  2  (A:,  q)  — 


8 


Ni 


p  n  k  +  1 


k 

x  (n,  p,  q)  —  £  afu)x  (n  —j,  p,  q) 
7“1 


1  + 


k 

z 

7-1 


(a) 


(2.6) 


for  u  —  1,  2,...  and  o-q  —  0.  Let  be  the  column  vector  whose  components  are  the 
approximate  values  of  the  coefficient,  a ,  ’s,  obtained  on  the  u-th  iteration  by  solving  the 
system  of  linear  equations 

(M(k,  q)—fl  <r2_j  (k,  q) /)  a(u)  —  X(k,q),  (2.7) 


where  I  is  an  identity  matrix.  In  summary,  we  choose  positive  integers  k  and  q  and  use 
the  tabular  data  x  (n  )  to  generate  the  augmented  matrix  (M(A:,  q),  X(k,  q))  correspond¬ 
ing  to  the  overdetermined  system  (2.5).  Then  by  starting  with  o-J  =  0  we  obtain  the  first 
approximations,  af1'*' s,  by  solving  (2.7).  These  values  are  substituted  in  (2.6)  to  obtain 
o-2  which  in  turn  is  used  in  (2.7)  to  obtain  the  second  approximation  vector  a(2\  and  so 
on.  The  iteration  is  continued  until  |  —  aj-u^  |  <  8  —  10-6  or  u  —  20.  Together 

with  a  solution  of  (2.7)  the  determinant  of  this  system  is  also  computed.  The  absolute 
value  of  the  determinant  of  the  last  iteration  is  denoted  by  D  (k,q ).  Thus,  at  the  end  of 
the  iteration  we  have 


computed  by  (2.6), 


and 


o-2(k  ,q  ) 

0*1  ( k,q\  a 2  (k ,q ),...,  ak  (k,q)\ 
D  ( k,q ). 


The  corresponding  autoregressive  model  is 

x  ( nj>,q )-  £  aj  ( k,q)x  ( n—j,p,q ). 
7-1 

The  eigenvalues  of  the  autoregressive  model,  namely 


(2.8) 

(2.9) 

(2.10) 

(2.11) 
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X  l(k,qX  X  2(k,q\...,  X  k(k,q ), 


(2.12) 


are  the  roots  of  the  equation 

\k  (k,q)—  £  dj  (k,q) \k~J  (£,<?).  (2.13) 

J  -  1 

It  is  very  unlikely  that  the  estimates  of  the  aj  ’s  obtained  from  the  tabular  data  will  yield 
multiple  roots  of  (2.13).  Therefore,  we  assume  here  that  all  the  roots  are  distinct,  since  all 
the  cases  of  synthetic  data  examined  in  this  report  and  empirical  data  of  the  previous 
report  2  produced  no  pairs  of  nearly  equal  roots.  We  note  that  exact  representation  of 
y  (t)  =4' (1  4- 1  +t2  + t3)  used  below  to  generate  one  set  of  synthetic  data  does  yield  four 
equal  roots.  However,  an  approximation  of  the  values  of  this  function  with  an  added 
noise  does  not  produce  an  autoregressive  model  with  equal  or  nearly  equal  eigenvalues. 

As  stated  above,  an  approximating  function  is  an  element  of  the  algebra  A.  Hence, 
with  simple  eigenvalues  we  take 

x  (n,p,q)~  £  CjXf  (2.14) 

for  all  data  points  with  the  same  value  of  n .  The  coefficients,  Cj ’s,  are  selected  to  minim¬ 
ize  the  RMSE  of  the  resulting  approximation  over  a  span,  say,  from  n—K  to  n+K.  We 
select  K^k.  Thus,  the  coefficients,  Cj  ’s,  are  obtained  by  solving  the  normal  equations  of 
the  overdetermined  system  (2.14),  i.e.,  by  solving  the  following: 

£  i^X/^X"4^-  £  x(n+u,p,q)X^,i  -  1,2,...,A:.  (2.15) 

u  —  —K  7—1  u  —  —K 

k 

Let  Cl  —  (clt  c2,  •  •  ■  .ck),  Sjj  —  £  X/'X",  S  be  the  matrix  (s^),  and  Z  be  the 

ku--k 

column  vector  with  components  £  Xj^x  (n+u,p,q  ),  J  —  1,2, ...,k.  Then  (2.15)  can 

u  —  —k 

be  rewritten  in  equivalent  matrix  form:  S  C  k  =  Z.  Consequently, 

Ck-S~lZ  (2.16) 

Let  Y„  be  the  column  vector  with  components  x  ( n  +  u,p,q)\f,  j  —  1,2 ,...,k.  Then 
we  get  from  (2.16): 

Ck-  £  S~lYu.  (2.17) 

u  —  —k 

By  writing  \T  =(  Xf,  X£,  .  .  .  ,X£}  and  by  substituting  (2.17)  in  (2.14)  we  get  the 
smoothed  value  of  the  function: 

xk{n,p,q)~  £  A  tS~'Yu.  (2.18) 

u  —  — k 


10 


In  view  of  (2.3)  the  relation  (2.14)  can  be  written  as  follows: 

St  *(*)=  £  Cj  A  q  ,  or 
7-1 

4(r+pA)-  T.  ejxf~E~  •  (219) 

7-i 

Differentiation  of  (2.19)  5  times  yields 

d’*U  +p  A) - 1—  V  c.X/^log’X,  .  (2.20) 

We  replace  f  in  (2.20)  by  t  —p  A,  and  then  we  let  f  —  r  A  with  r  as  in  (2.3).  This  in 
turn  leads  to 

(njfrf)'-  log**,.  (2.21) 

In  view  of  (2.17)  the  estimate  of  the  s-f/t  derivative  in  (2.21)  can  be  written  as  fol¬ 
lows: 

£  A/S-%  (2.22) 

(q&r  u  mm-k 

where  \s  —  ds  A  /  du  s  is  a  column  vector  with  components  AjMog *\jJ  —  1,2,..., A: 
This  completes  the  description  of  the  basic  relations  of  the  algorithm. 


3.  RELATIVE  ACCURACY  OF  APPROXIMATIONS 


The  algorithm  of  the  preceding  section  yields  several  approximating  functions 
corresponding  to  various  choices  of  k  and  q .  According  to  the  choices  in  an  earlier 


report2  there  are  3 


min(-jy,39) 


approximating  functions  for  a  data  set  of  N  points. 


Here  the  square  brackets  denote  an  integer  not  exceeding  the  expression  in  the  brackets. 
We  discuss  now  a  heuristic  rule  to  assign  a  weight  to  each  approximation. 


As  described  in  Section  2  of  ND  all  the  approximations  corresponding  to  autoregressive 
models  with  eigenvalues  containing  a  real  negative  part  are  assigned  weight  zero,  since  in 
this  case  an  approximating  function  contains  a  periodic  term  with  the  period  so  short  that 
less  than  four  data  points  are  contained  in  the  period;  i.e.,  the  data  are  inadequate  to  deter¬ 
mine  this  term  with  sufficient  accuracy.  We  illustrate  this  by  an  approximation  to  x  (r)  = 
erf  (r )  +  c,  where  e  is  a  normally  distributed  random  variable  with  standard  deviation  tr. 
We  generate  data  by  evaluating  x  (f )  in  the  interval  [0,2]  with  the  step  size  A  =  0.005 
and  by  adding  pseudorandom  numbers  e.  Autoregressive  models  computed  for  the  values 
of  this  function  plus  random  error  with  o- —  0.00119  corresponding  to  eigenvalues  with 
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positive  real  parts  allow  a  maximum  frequency  to  max  —  0.45ir.  If  only  the  terms  with  the 
i  7 r 

frequencies  to ,  —  —  satisfying  the  condition  to  l_l<tomllx^tol  are  included  in  a  trun¬ 
cated  Fourier  series  we  obtain  an  approximation  with  RMSE’s  equal  to  0.087,  0.420,  and 
1.10  for  the  function,  first,  and  second  derivatives,  respectively.  If  the  terms  of  the  next 
higher  frequency  are  added  to  this  approximation,  the  RMSE’s  are  0.068,  0.537  and  2.41. 
Thus,  we  have  a  slight  improvement  in  the  v'-es  of  the  function,  as  expected,  but  more 
than  twice  as  big  an  error  in  the  second  derivative.  This  example  shows  that  rejection  of 
the  autoregressive  models  with  negative  parts  in  their  eigenvalues  is  appropriate. 

We  made  similar  tests  in  two  additional  cases,  namely,  with  the  data  obtained  by 
evaluating  the  Bessel  function  J0(t )  of  the  first  kind  (500  points  in  the  interval  [1,6]  with 
one-sigma  error  of  0.00344)  and  by  evaluating  sin27r  t  +  0.1sinl0  n  t  (250  points  in  the 
interval  [0,  1]  with  one-sigma  error  of  0.022). 

In  the  first  case,  if  only  the  terms  with  the  frequencies  toj^toj  where 
a>,_i  <6>max^  to  j  are  included  in  the  Fourier  expansion  the  corresponding  errors  are 
0.121,  0.255,  and  0.334.  If  the  terms  with  frequency  to  /+1  are  added  the  errors  are  0.090, 
0.296,  and  0.588.  Again  we  have  a  slight  improvement  in  the  values  of  the  function,  but 
derivatives  become  less  accurate. 

In  the  second  case  the  terms  with  to  ,J  =  1,2, 3,4,5  are  allowed  in  the  Fourier  expan¬ 
sion,  according  to  the  same  criterion.  The  resulting  RMSE’s  in  the  values  of  the  function, 
first,  and  second  derivatives  are  0.0038,  0.074,  and  1.73.  If  the  terms  with  the  frequency 
to  5  are  dropped,  the  RMSE’s  become  0.07013,  2.23,  and  69.6,  respectively.  If  the  terms 
with  tojJ  =  1,2, 3, 4, 5, 6  are  included  the  errors  are  0.01169,  0.209,  and  6.09.  Thus,  here 
again  the  selected  criterion  for  cut-off  frequency  in  the  representation  of  the  data  yields 
optimal  results. 

Similarly,  comparison  of  approximations  computed  by  the  method  of  the  preceding  sec¬ 
tion  corresponding  to  autoregressive  models  with  negative  real  parts  and  those  correspond¬ 
ing  to  positive  real  parts  shows  that  rejection  of  approximations  with  negative  real  parts  in 
their  eigenvalues  is  appropriate. 

For  instance,  for  x  (f )  **  sin27rr  (with  noise)  evaluated  in  [0,  1  ]  at  a  step  of  0.004  the 
constraints  imposed  on  k  and  q  allow  a  total  of  48  models.  Of  these,  24  contain  eigen¬ 
values  with  negative  real  parts.  The  most  accurate  approximation  among  these  yields 
RMSE’s  of  0.069  and  2.17  in  the  first  and  second  derivatives,  while  the  best  approximation 
with  positive  real  parts  in  eigenvalues  yields  0.0397  and  0.0920,  respectively. 

The  remarks  in  Section  2  of  ND  and  the  examples  just  presented  provide  heuristic  and 
empirical  justification  for  rejecting  (i.e.,  assigning  weights  of  zero  to)  the  approximations 
corresponding  to  autoregressive  models  with  negative  real  parts  in  their  eigenvalues. 

2 

The  weights  assigned  in  a  previous  report  to  the  approximations  with  positive  real 
parts  in  their  eigenvalues  produced  relatively  good  results.  Here  we  examine  a  modified 
form  of  the  weights  that  yields  somewhat  better  results.  We  assume  as  before  that  the 
weight  increases  as  the  RMSE  cr(k,q )  of  the  solution  of  the  overdetermined  system 
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(2.5)  decreases  and  that  it  also  increases  as  the  corresponding  normal  system  becomes  less 
sensitive  to  perturbations  of  the  data,  i.e.,as  the  absolute  value  of  det  M  (k,q)  in  (2.10) 
increases.  We  select  a  simple  function  satisfying  these  conditions  and  invariant  with 
respect  to  rescaling  of  the  dependent  variable  x ;  i.e., we  choose  the  weight  w  ( k,q )  given 
by 


w  (k,q) 


D  ( k,q ) 
o^ik.q) 


(3.1) 


Comparison  of  accuracy  of  various  models  provides  empirical  support  for  this  choice  of 
their  relative  weights.  Eight  sets  of  synthetic  data  yield  a  total  of  1%  approximating  func¬ 
tions  for  various  values  of  k  and  q  with  non-negative  real  parts  of  eigenvalues.  The 
number  of  the  models  of  this  kind  for  individual  data  sets  range  between  16  and  37. 
Thus,  we  have  196x2—392  approximations  of  the  first  and  second  derivatives.  It  turns  out 
that  the  three  highest  weights  computed  by  (3.1)  for  each  set  of  data  differ  very  little,  and 
the  values  of  the  weights  drop  rapidly  for  the  remaining  models.  For  this  reason  we  con¬ 
sider  the  three  models  with  the  highest  weights  for  each  set  of  data.  We  have  altogether 
48  approximations  of  the  first  and  second  derivatives  with  the  weights  among  the  three 
highest.  Among  these  there  are  32  approximations  that  produce  smaller  errors  in  the 
respective  derivatives  than  the  remaining  360  approximations.  Nine  out  of  sixteen  approx¬ 
imations  with  the  highest  weights  in  their  respective  sets  actually  are  best  for  the  first  or 
second  derivative.  The  remaining  seven  models  with  highest  values  of  w  (k,q)  do  not 
yield  best  approximations.  However,  they  are  not  much  inferior  to  the  best  models  of  the 
autoregressive  type.  In  fact,  the  average  reduction  of  the  RMSE  of  derivatives  between 
the  models  with  the  highest  w  (k,q)  and  the  best  ones  by  comparison  with  exact  values 
of  derivatives  is  only  32%  for  these  seven  cases.  The  results  of  the  next  section  show  that 
the  RMSE  of  derivatives  computed  by  a  spline  approximation  is  higher  in  one  case 
(second  derivative  of  ef)  by  as  much  as  26700%  than  the  RMSE  obtained  by  the  current 
method.  Even  if  this  extreme  case  is  ignored,  the  RMSE  of  the  approximations  computed 
by  the  spline  procedure  is  474%  higher  than  that  of  the  current  method.  In  view  of  this  a 
32%  increase  in  those  cases  in  which  the  criterion  (3.1)  fails  to  select  the  best  model  is  not 
significant.  Consequently,  we  accept  this  criterion  on  the  basis  of  empirical  results. 

The  smoothed  functional  values  and  estimates  of  derivatives  for  each  ( k  ,q  )  pair  are 
given  by  (2.18)  and  (2.22),  respectively,  with  p  —  0,1,2,  •••  ,<?— 1  and 
n  —  fc+1,  fc+2,  •  •  ■  ,Np—k ,  where  Np  is  the  same  as  in  (2.6).  This  produces  no 
smoothed  values  for  t  —  A, 2  A  ,  •  •  •  ,  kq  A  ,  and  for 

t  —  (N  —kq  4-  1)A,(N  —  kq  4-  2)A,  •  ■  •  ,N  A.  Denote  the  values  of  k  and  q 
corresponding  to  the  three  models  with  highest  values  of  w  (k,q)  by  (/cy  qt\i  —  1,2,3  and 
let  n0  —  max  ( ktq We  obtain  the  smoothed  values  and  approximation  of  derivatives 
based  on  all  three  models  with  highest  weights  only  at  the  points  t  —  (n0  +1)A, 
(n0  —  2)  A,  •  •  ■  ,(iV  —  n0  )A.  These  points  belong  to  what  we  call  here  the  central  inter¬ 
val  of  the  data.  The  smoothed  functional  values  %  («)  and  approximate  derivatives 
Xis\n)  at  the  points  of  the  central  interval  are  defined  by  the  following: 
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(3.2) 


2  (r )  ■=£*>  (k ,q )  Xk  (n  j),q)  / £w  (fc,g) 

and 

x(s)(r)=2>  a,9)x(s)  {n,p,q)  (k,q)  (3.3) 

k,q  k,q 

with  r  —  p  4-  qn. 

Instead  of  estimating  derivatives  at  the  data  points  not  included  in  the  central  interval 
with  the  aid  of  available  models  of  lesser  accuracy  (such  as  with  k  —  q  —  1)  we  adopt  the 
filtering  procedure  described  in  the  next  section. 


4.  DERIVATIVES  IN  THE  INITIAL  AND  FINAL  SEGMENTS  OF  THE  DATA 

Let  k0  —  max  (fcj,  k2,  k  3),  where  k{  ,  /  —  1,2,3  are  the  parameters  of  the  three  models 
with  the  highest  weights  w  (/:,  ,  qt ),  and  let  q0  be  the  value  of  q, ,  such  that  ki  —  k0 . 
Equations  (3.2)  and  (3.3)  yield  k  (r)  and  i (r)  "for  r  *■=  (n0  +  1 )  A, 
(  n0  4-  2)  A,  ■  •  •  ,(/r  —  n0)  A.  We  assume  digital  filtering  models: 

*'C/)-  Z  di*  0  ~  iq0 )  +  Z  /,  x'  (j  —  iq0)  (4.1) 

/  -  o  /  -  1 

for  each  pair  of  non-negative  integers  p0  and  r0  satisfying  the  conditions  pa  ,  rn  . 


If  j  is  such  that  n0  +  q0  max  (p0  ,  r0)  ^  N  —  n0  then  the  values 

St'  (j)  ,  St  (J  —  iq0)  and  St  '  (J  —  iqQ )  belong  to  the  central  interval  of  the  data.  Hence, for 
these  values  of  j  (4.1)  constitutes  an  overdetermined  system  with  unknowns  and  /,  . 
We  solve  this  system  by  the  least  squares  method  and  obtain  the  estimates  d,  ( r0  ,  p0  ) 
and  fj  (r0  ,  p0 )  ,  the  RMSE  o-  (p0  ,r0 )  of  the  solution,  and  the  absolute  value 
E  (p0  ,r0)  of  the  determinant  of  the  corresponding  normal  equations.  In  analogy  to  the 
preceding  section  we  assign  a  relative  weight  W  ( p0  ,rc)  to  each  of  the  digital  filters  (4.1) 
as  follows: 


W  (. P0S0 ) 


E  ( P0J0 ) 


o-0Vro> 


ft  +1 


(4.2) 


Let  (Pi,r  j),  (p2  ^2)  and  0?3  ,r3)  be  the  pairs  of  (p0r0 )  corresponding  to  the  three  largest 
values  of  W  (p0  ,ra  )  .  For  s  —  1,2,3  let 

Xj  ( ps  ,rs)  0)  , 

if  n0  +  1  <  j  ^.N  —  n0  .  Let 

Xj  ips  ,rs )  -  x  0)  , 


ify>iV  —  n0  .  Let 
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*/  ips  ,r5)  —  %'  0)  » 


if  n0  +  1^  j  <  N  —  n0  .  Let 

*  ' ,(ps,rs )  -  £  </,  (ps  ,rs)x(J  -  iq0  )  4-  £  ft  (ps  ,rs)2‘  (J  -  iq0)  , 

/ -o  / - 1 

if  N  >N  —  n0  .  Then  the  estimates  of  the  first  derivative  in  the  final  interval  of  the 
data  are  given  by 

*'(/)-  £  W  (ps ,  rs)%j  (ps,rs)/  Z  W  (ps  ,rs)  (4.3) 

S  —  1  5—1 

j  -  TV  -  n0  +  1  ,N  -  /I0  +  2 . N. 

The  estimates  of  the  second  derivative  are  obtained  in  the  same  manner  by  replacing 
jt'ij)  by  *"0‘)  for  j  =  n0  +  1,  n0  +  2  ,  •  •  •  ,  N  —  n0  in  (4.1)  and  then  by  solving  the 
resulting  overdetermined  system,  i.e.,  by  determining  the  corresponding  linear  digital 
filters.  These  are  used  to  obtain  the  estimates  of  the  second  derivative  analogous  to  (4.3). 

The  estimates  of  the  derivatives  in  the  initial  segment  of  the  data  are  obtained  by 
renumbering  the  data  points  as  follows:  £  0 )  —  x  (N  —  j  +  D  and  then  by  applying  the 
procedure  just  described  to  £  0  )  • 

5.  NUMERICAL  EXAMPLES 

We  compute  approximations  of  the  first  and  second  derivative  by  the  method  described 
on  the  preceding  pages  for  ten  sets  of  synthetic  data.  The  ten  cases  numbered  in  Column 
1  of  Table  1  are  obtained  by  evaluating  the  functions  x  (t)  listed  in  Column  2  in  the 
intervals  given  in  Column  3  at  the  step  size  A  shown  in  Column  4  and  then  by  adding  to 
each  value  a  pseudo-random  error  normally  distributed  with  zero  mean  and  standard  devi¬ 
ation  o-  as  contained  in  Column  5.  In  each  case  the  standard  deviation  cr  is  equal  to  the 
average  of  the  absolute  values  of  the  change  of  the  corresponding  function  as  its  argument 
changes  by  A  .  With  random  errors  of  this  size  derivatives  cannot  be  estimated  by  divided 
differences. 

In  Case  9  cr  —  0  ,  i.e.,  here  we  have  functional  values  exact  up  to  15  decimal  digits, 
which  is  single  precision  for  the  computer  employed. 

In  Case  5  the  function  x  (f )  —  J0  (f )  is  the  Bessel  function  of  the  first  kind.  In  Case  7 

(Witch  of  Agnesi)  x  (f)  —  2  cos2  (tan_1y  )  and  in  Case  8  we  have  a  rational  function 

x  (f )  —  36  t  /  (f2  +  9)  .  These  three  functions  as  well  as  the  error  function  are  included 
to  test  how  well  their  derivatives  can  be  approximated  when  the  functions  do  not  belong 
to  the  Algebra  A  defined  in  Section  2. 
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TABLE  1.  SYNTHETTC  DATA 


1 

2 

3 

4 

5 

Cue 

*  (l) 

Tntervtl 

A 

<r 

1 

sin  2  » t 

[0,1] 

.004 

.001 

2 

sin  2  it  t  +0.1  sin  10  *  I 

[0,1] 

.004 

.022 

3 

<  ' 

[0,5] 

.01 

.295 

4 

4'  (1  +1  +  <2+l3) 

[0v2] 

.004 

.478 

5 

/.(») 

[1.6] 

.01 

.00344 

6 

error  function 

[0,2] 

.005 

.00119 

7 

Witch  of  Agnesi 

[0^] 

.005 

.0025 

8 

Newton  Serpentine 

[0,2] 

.005 

.0139 

9 

2  r3  —  9f2  +  12  r 

[0,31 

.005 

0 

10 

2r3— 9rJ  +  l2r 

[0,3] 

.005 

.01833 
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The  results  are  summarized  in  Table  2.  Here  Column  1  is  the  same  as  in  Table  1.  To 
each  row  of  Column  1  there  correspond  two  rows  of  the  remaining  columns.  One  of  these 
contains  the  errors  of  the  first  derivative  x'(t )  and  the  other  that  of  the  second  derivative 
as  indicated  in  Column  2.  Columns  3  and  4  give  RMSE’s  of  the  respective  esti¬ 
mates  of  derivatives  expressed  in  percent  of  the  RMS  of  the  actual  values  at  the  data 
points.  Column  3  (B.S.)  shows  the  errors  of  the  estimates  obtained  by  fitting  the  noisy 
data  by  B-splines.  Column  4  contains  the  errors  (in  percent)  of  the  estimates  computed  by 
the  current  method  (C.M.),  i.e.,  by  the  method  of  this  report.  By  comparing  these  two 
columns  we  see  that  the  method  described  here  gives  much  better  estimates  of  derivatives 
in  all  but  three  cases  out  of  20.  Two  of  these  cases  are  for  exact  values  of  the  cubic  poly¬ 
nomial.  The  first  and  the  second  derivatives  computed  by  cubic  B-spline  approximation 
agree  up  to  eight  digits  with  the  exact  values  in  this  case.  The  coefficients  of  the  autore¬ 
gressive  model  with  highest  weight  agree  with  the  exact  coefficients  up  to  seven  digits. 
The  resulting  perturbation  of  one  or  two  units  in  the  eighth  digit  yields  eigenvalues  accu¬ 
rate  up  to  only  three  decimal  digits.  The  derivatives  computed  with  these  approximations 
on  the  average  differ  from  the  exact  values  in  the  fourth  decimal  digit,  as  shown  in 
Column  4.  The  last  two  rows  show  that  the  method  of  this  report  produces  considerably 
better  results  than  the  B-spline  approximation  when  the  values  of  this  same  polynomial 
contain  random  errors. 

In  Case  2  the  error  for  the  first  derivative  obtained  by  the  method  of  this  report  is 
almost  twice  as  big  as  that  obtained  by  the  spline  approximation.  However,  the  error  of 
the  second  derivative  computed  by  the  current  method  is  only  half  the  error  obtained  by 
the  spline. 

Of  the  remaining  sixteen  derivatives  (aside  from  Cases  2  and  9)  five  approximations  by 
the  current  method  are  better  by  an  order  of  magnitude  or  more  than  that  by  the  spline; 
six  approximations  have  errors  four  or  more  times  smaller  than  the  errors  by  the  spline 
method,  and  the  remaining  five  are  at  least  twice  as  good  as  those  by  the  spline  procedure. 
The  worst  error  by  the  spline  approximation  as  compared  to  the  approximation  by  the 
current  method  is  that  of  the  second  derivative  of  the  exponential  function.  Even  if  this 
worst  case  is  ignored,  the  error  by  the  spline  approximation  is  on  the  average  474%  of  the 
error  by  the  method  of  this  report. 

Compatibility  of  derivatives  computed  by  our  method  with  the  data  was  further  tested 
by  comparing  exact  functional  values  with  an  approximate  first  integral  of  the  first  deriva¬ 
tive  and  the  second  integral  of  the  second  derivative.  Approximations  of  the  integrals 
were  obtained  by  trapezoidal  rule.  The  RMSE  of  the  integral  of  the  first  derivative  com¬ 
puted  for  various  segments  of  the  central  portion  of  the  data  varies  from  0.075  to  1.88  per¬ 
cent  of  the  RMS  of  the  functional  values.  The  smallest  error  is  for  Witch  of  Agnesi  (Case 
7)  and  the  largest  for  Case  4.  The  values  of  the  derivatives  in  the  initial  and  final  data 
segments,  however,  are  less  accurate.  Consequently,  the  RMSE’s  of  the  integrals  of  the 
first  derivative  computed  for  various  segments  of  the  complete  data  set  vary  from  0.10 
(Case  3)  to  3.8  (Case  5)  percent.  Similarly,  the  RMSE’s  of  the  second  integral  of  the 
second  derivative  in  the  oentral  portion  of  the  data  vary  from  0.26  (Case  3)  to  17.5  (Case 
4)  percent  and  in  the  complete  set  of  data  between  0.25  (Case  3)  and  43.5  (Case  2)  per¬ 
cent. 
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TABLE  2.  APPROXIMATION  ERRORS  IN  PERCENT 


1 

2 

3 

4 

Case 

Deriv 

B.S. 

C.M. 

x' 

2.5 

.61 

1 

x" 

20 

1.1 

x' 

3.3 

6.2 

2 

X" 

23 

11.9 

x' 

3.3 

.22 

3 

x" 

67 

.25 

x' 

3.2 

1.6 

4 

x" 

67 

8.6 

x' 

4.7 

2.2 

5 

x" 

140 

15.5 

x' 

4.5 

1.7 

6 

x" 

253 

37 

x' 

3.8 

1.9 

7 

x" 

191 

36 

x' 

4.7 

1.1 

8 

X" 

347 

28 

x' 

0 

.042 

9 

x" 

0 

.071 

x' 

4.29 

1.55 

10 

x" 

88 

3.55 
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We  note  further  that  our  method  is  not  too  time-consuming  in  spite  of  the  fact  that  it 
determines  several  approximations  and  the  corresponding  computer  program  contains 
several  loops.  Computing  time  on  the  CDC  CYBER  76  computer  of  BRL,  including  gen¬ 
eration  of  synthetic  data,  in  none  of  the  ten  cases  exceeded  8.8  sec.  The  average  comput¬ 
ing  time  per  case  was  4.5  sec. 


6.  EXTRAPOLATION/FORECASTING 

A  variant  of  the  new  method  described  in  this  and  earlier  related  reports  can  be  used 
for  extrapolation  out  of  the  region  in  which  the  data  was  sampled.  As  with  any  forecasting 
procedure  it  must  be  assumed  that  the  underlying  trend  shown  by  the  sample  continues 
into  the  region  considered.  If  conditions  change  that  trend  significantly,  the  forecast  will 
naturally  be  much  less  reliable. 

If  equation  (2.14)  remains  valid,  functional  values  s  ■  q  ■  A  t  later  in  time  are  obtained 
by  replacing  n  by  n  4-  5  ,  while  the  eigenvalues  A  remain  the  same.  Other  equations  are 
modified  accordingly. 

In  the  central  portion  of  the  sample  a  span  of  2  k  +  1  points,  q  ■  A  t  apart,  was  used  to 
find  values  at  a  point  t  offset  zero  from  the  center  of  the  span.  In  the  right-hand  end  por¬ 
tion  the  offset  was  changed  to  k  ■  q  ■  A  t  so  as  to  use  available  data  points,  to  the  left  of 
and  including  t ,  in  the  evaluation.  If  the  offset  is  increased  further,  to  v  ■  q  •  A  t  for 
some  v>k,  the  span  can  remain  within  the  data  sample  even  though  t  is  outside.  This 
provides  the  basis  of  forecasting. 

The  correspondingly  modified  algorithm  was  used  for  extrapolation  of  several  of  the 
standard  test  functions.  In  cases  where  the  A ’s  closely  approximate  those  of  the  exact 
(non-noisy)  function,  extrapolated  values  well  outside  the  sample  were  obtained  with 
acceptable  accuracy.  In  cases  of  greater  noise,  poorer  approximation  of  the  true  A ’s  is 
typical,  and  extrapolation  is  therefore  less  satisfactory. 

For  x  =  sin  2  rr  t  +  e  ,  201  points  on  [0,1],  prediction  of  the  extrema  and  zeroes  on 
[1,2]  was  undertaken.  For  a  noise  level  (  cr  (  e  »of  .001  the  zeroes  at  t  =  1.5  and  2.0 
were  predicted  to  be  at  t  dr  1  •  A  t  with  error  in  St  (t )  less  than  .001  in  those  intervals. 
Similarly,  the  extrema  at  1.25  and  1.75  were  predicted  at  t  ±  1  •  A  t  with  |  St  |  values  on 
(0.999,  1.000).  The  corresponding  first  derivatives  computed  by  the  linear  filtering 
method  of  Section  4,  were  found  with  RMS  error  within  0.2  percent  of  the  overall  RMS 
(x). 

Corresponding  results  with  greater  noise  (  cr  =.010)  also  located  the  extrema  and 
zeroes  within  1-  A  t,  with  |£|  in  error  by  .010  (RMS  error)  over  the  entire  interval  [1,  2] 
The  RMS  error  for  x'  in  the  extrapolated  region  was  3.0%  of  RMS  (x'  ). 

For  the  more  oscillatory  function  x  =  sin2  n  t  +  0.1  sin  10  rr  t  4-  e  ,  with 
tr  (  e  )  =.001,  the  extrema  at  t  =  1.25  and  1.75  and  the  zeroes  at  t  =1.5  and  2.0  were 
predicted  to  within  1  •  A  t,  and  the  values  near  the  extrema  were  correct  within  .002.  The 
corresponding  first  derivatives  St'  near  the  zero  at  t  =1.5  were  in  error  by  less  than  3%  of 
true  value,  degrading  somewhat  thereafter.  The  corresponding  case  with  cr  =  .010 
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extrapolated  j?  well  (RMSE  of  .02  on  [1,  2])  but  x'  began  to  degrade  somewhat  earlier  in 
the  extrapolation  region.  Note  that  the  increased  input  noise  is  a  significant  fraction  of  the 
amplitude  of  the  second  term  of  this  double  sine  function. 

Extrapolation  should  always  be  used  with  caution,  but  there  are  many  practical  cases 
where  it  may  be  necessary.  Economic  forecasting,  demographic  projections,  and  prediction 
of  future  locations  of  moving  targets  are  examples.  The  reverse  process  of  extrapolating 
into  the  past  is  also  of  interest,  e.g.,  likely  earlier  values  of  a  variable  with  only  recent 
measurements  available.  This  can  be  achieved  by  storing  the  data  values  in  reverse  and 
changing  the  sign  of  first,  third,  and  other  odd  derivatives  of  interest.  The  current  gen¬ 
eralized  algorithm  thus  can  be  used  in  either  context. 
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APPENDIX 


DNDLIB  -  A  USER  LIBRARY  FOR  DIFFERENTIATING  NOISY  DATA 

There  now  exists  on  MFZ,  the  BRL  CDC  Cyber  7600  computer,  a  user  library  called 
DNDLIB,  created  by  G.  Francis  in  September  1981  and  based  on  the  computer  programs 
and  theoretical  algorithms  discussed  in  three  BRL  publications  by  C.  Masaitis  and  G. 
Francis,  this  report  and  its  earlier  companions. 

The  library  was  designed  for  easy  access  by  other  users  of  the  BRL  computer.  In  gen¬ 
eral  only  one  of  the  14  subroutines  included  need  be  called  by  a  user  program,  as  the  oth¬ 
ers  are  called  as  needed  automatically.  That  particular  subroutine  is  named  DND  and 
allows  a  certain  amount  of  flexibility  on  the  part  of  the  user.  He  must  set  up  an  array 
(vector)  of  equally  spaced  data  points  (usually  experimental  values  containing  noise,  i.e., 
measurement  errors  and  the  like)  and  specify  the  number  of  such  points,  say  N.  In  addi¬ 
tion  he  must  select  the  maximum  order  of  derivatives  of  interest  (1  to  3)  and  provide 
additional  arrays  for  smoothed  results  and  each  order  of  derivative.  The  current  upper 
bound  on  N  is  1000,  and  N  should  exceed  30. 

A  sample  call  (in  Fortran)  is  as  follows: 

CALL  DND  (XN,  N,  DT,  MXOD,  XS,  Dl,  D2,  D3) 

where,  say,  N  =  250,  XN  is  an  array  of  at  least  N  values  (points  1  to  N  to  be  processed), 
DT  is  the  spacing  of  the  independent  variable,  MXOD  is  the  maximum  order  of  derivative 
wanted  (1  to  3),  and  arrays  XS  (of  length  N  or  more)  and  Dl  (similarly)  are  to  receive 
the  smoothed  values  and  the  first  derivatives  at  corresponding  points.  If  MXOD  is  2  or 
more,  an  array  D2  must  be  provided;  likewise  if  MXOD  is  3,  an  array  D3  is  needed.  The 
labels  D2  and  D3  may  be  omitted  from  the  call  if  MXOD  is  1 . 

For  each  point  i,  i=*l  to  N,  there  will  be  a  smoothed  value  at  XS(i)  and  a  first  deriva¬ 
tive  at  Dl(i).  If  MXOD  is  2  a  numerical  second  derivative  will  be  at  D2(i),  and  if  MXOD 
is  3  then  D3(i)  will  contain  the  third  derivative  corresponding  to  point  i.  (If  still  higher 
derivatives  are  desired,  the  process  can  be  repeated  by  copying  D3,  say,  into  XM  and  mak¬ 
ing  a  new  call.  The  accuracy  of  derivatives  grows  worse,  however,  as  the  order  increases, 
particularly  if  the  original  data  contains  considerable  error.) 

The  method  of  DNDLIB  uses  a  few  subroutines  from  the  IMSL  library,  so  that  library 
must  be  made  available,  too.  This  is  done  on  MFZ  by  means  of  the  following  CDC  ’con¬ 
trol  cards’  or  suitable  replacements: 

ATTACH(IMSL) 

ATTACH  (DNDLIB,  ID=PUBLIC) 

LIBRARY(IMSL,DNDLIB) 

In  addition  the  standard  BRL  subroutines  known  as  FNEQS  and  MATINV  are  used  for 
matrix  setup  and  inversion.  These  are  available  on  MFZ  with  no  user  action  required. 
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The  size  of  the  full  set  of  subroutines  is  approximately  20K  words  of  SCM  (central 
memory)  and  12K  of  LCM  plus  the  user-supplied  arrays  (of  3N  to  5N  words,  as  men¬ 
tioned  earlier).  This  includes  all  work  space  for  N  up  to  1000  and  MXOD  up  to  3,  mostly 
in  LCM. 

The  time  required  for  calculations  is  highly  dependent  on  the  number  of  data  points 
and  also  on  the  order  of  derivatives  wanted  as  well  as  the  degree  of  complexity  of  the 
underlying  function.  For  a  sample  noisy  sinusoid  with  derivatives  wanted  through  order  2 
the  following  times  were  found  (in  cpu  seconds  on  MFZ): 

N  -  100  250  500  1000 
T-0.5  1.3  4.4  11.1 

As  indicated  above  results  are  found  at  all  N  points,  not  just  a  few. 

If  extrapolation  is  required,  an  alternate  call  of  much  the  same  form  is  used,  with  one 
additional  parameter,  denoted  NP,  and  the  subroutine  name  is  DNDE  rather  than  DND: 

CALL  DNDE  (XN,  N,  DT,  MXOD,  NP,  XS,  Dl,  D2,  D3) 

A  request  for  extrapolation  forward  is  entered  by  specifying  NP  greater  than  N,  but  under 
1000.  This  feature  should  be  used  with  caution.  Naturally,  all  arrays  for  results  should  be 
of  at  least  NP  cells. 

It  is  hoped  that  this  library  will  prove  of  benefit  to  users,  who  are  encouraged  to  apply 
it  to  their  own  noisy  data,  especially  when  the  assumptions  of  Section  2  are  thought  to  be 
satisfied. 
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